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We address dissipation effects on the non-equilibrium quantum dynamics of an ensemble of spins- 
1/2 coupled via an Ising interaction. Dissipation is modeled by a (ohmic) bath of harmonic oscillators 
at zero temperature and correspond either to the sound modes of a one-dimensional Bose-Einstein 
(quasi-)condensate or to the zero-point fluctuations of a long transmission line. We consider the 
dimer comprising two spins and the quantum Ising chain with long-range interactions, and develop a 
(mathematically and numerically) exact stochastic approach to address non-equilibrium protocols in 
the presence of an environment. For the two spin case, we hrst investigate the dissipative quantum 
phase transition induced by the environment through quantum quenches, and study the effect of 
the environment on the synchronization properties. Then, we address Landau-Zener-Stueckelberg- 
Majorana protocols for two spins, and for the spin array. In this latter case, we adopt a stochastic 
mean-held point of view and present a Kibble-Zurek type argument to account for interaction effects 
in the lattice. Such dissipative quantum spin arrays can be realized in ultra-cold atoms, trapped 
ions, mesoscopic systems, and are related to Kondo lattice models. 


I. INTRODUCTION 

Spin-boson models play a major role in various 
branches of physics, from condensed-matter physics, 
quantum optics, quantum dissipation, to quantum 
computatiorP®. A large collection of harmonic oscil¬ 
lators (bosons) can simulate dissipation, resulting in 
the celebrated Caldeira-Leggett modeP, giving rise to 
dissipation-induced quantum phase transitions observed 
in various contexts^. For example, a ohmic bosonic 
bath can be engineered through a long transmission line 
or a one-dimensional Luttinger liquicP^. An environ¬ 
ment can also affect the critical exponents associated 
with a phase transition such as the disordered-ordered 
transition in the quantum Ising chairP^^^. 

An impurity spin embedded in an environment also 
emerges as an effective model for strongly correlated 
quantum matter within dynamical mean-held theor^D^. 
The spin-boson model can be seen as a variant of the 
Caldeira-Leggett model where the quantum particle is 
a spin-1/2. The spin-boson model with an Ohmic bath 
exhibits a variety of rich phenomena such as a dissipative 
quantum phase transition separating an unpolarized 
(delocalized) and a polarized (localized) phase for the 
spin, as well as a coherent-incoherent crossover in the 
dynamical Rabi-type propertieJ^. This model is also 
intimately related to Ising models with long-range forces 
and to Kondo physic^^^E^. 

Several theoretical methods have been devised to study 
the dissipative spin dynamics for one spin in an ohmic 
bath such as the non-interacting blip approximatiod^, 
Quantum Monte Carlo (QMC) methods on the Keldysh 
contouiII^Jf^, and the time-dependent (TD) Numerical 
Renormalization Group (NRG) approaclP^^^, with 


recent progress done concerning the treatment of driving 
and quenche^^. Stochastic approaches have been 
developed both in the context of stochastic wavefunction 
approached or Stochastic Schrodinger Equation (SSE) 
methods on the density matrid^^H^. Stochastic Liouville 
equations were obtained for the density matrix in Refs. 

EHMl 

In this paper, we first consider a cluster of two spins 
in such a ohmic bosonic bath. The two spins are cou¬ 
pled through an Ising interaction. This model, which 
can be realized in ultra-cold atomd^^, reveals a dis¬ 
sipative quantum phase transition similar to the one- 
spin sit uation, but occurring at a smaller dissipation 
strengttP^®^^, which facilitates the application of nu¬ 
merical methods such as the SSE method in a large win¬ 
dow of the phase diagram. Using the Rabi-type dynam¬ 
ics of the spin system, we reproduce the phase diagram 
obtained using the NRG approacld and QMCd, show¬ 
ing the trustability of the SSE method. We also com¬ 
pute spin-spin correlations induced by the bath at long 
time, and compare our results with those obtained with 
a variational approacll^. We quantitatively address the 
occurrence of synchronization between the two spins, in 
relation with the spin-spin correlation function. Then, we 
investigate non-equilibrium quenched dynamics far in the 
polarized phase, which has not been discussed previously 
in the literature, and also Landau-Zener-Stueckelberg- 
MajoransP^J®! type interferometry for the dimer model. 
Next, we consider a quantum Ising spin chain with long- 
range forces allowing a mean-field treatment for the spin 
dynamics. The main aspect we explore concerns the 
extension of Kibble-Zurek type physicJ^^^^ induced by 
magnetic field gradients in time (Landau-Zener sweeps) 
in the case of an interacting spin ensemble subject to dis¬ 
sipation. Applying the stochastic procedure as well as a 
physical argument, we describe the interplay between in- 
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teractions between spins and dissipative effects from the 
bath on the well-known Landau-Zener formulsP^®^.We 
note that recent theoretical works have addressed similar 
questions regarding the effect of macroscopic dissi pation 
on the dynamical properties of quantum spin array^^^^. 
We also note recent experiments in ultra-cold atoms ad¬ 
dressing Kibble-Zurek type physicJ^. 


A. Model 

Hereafter, we focus on a system of M interacting spins 
(for the dimer M = 2 and for a spin array M +oo), 
which are coherently coupled to one common bath of har¬ 
monic oscillators: 


A ^ ^ 

^ E < + E E + h) f 

p=l p=l k 

“ ^E'^x+E^^^^bfe- (1) 

p^r k 

Here, with z/ = {x, y, zj are Pauli matrices related 
to the spatial site p and the Planck constant h is set 
to unity. At each site, the states \^z,p)^ corresponding 
to the two eigenstates of with eigenvalues ±1, define 
the two possible orientations of the spin. The long-range 
ferromagnetic Ising interaction can be enginee red in sys- 
tems of trapped ion^^^H^ and ultra-cold atomJ^^MIIlEni, 
It can also be the result of the Van der Waals interaction 
in Rydberg medisP^^^^. This model can also be seen as 
an example of Kondo lattices in one dimension through 
bosonizatiorP^ (for a review on Kondo lattices, see for 
example Ref. [64D. 


B. Bath effects 

The interaction with the bath plays an important 
role and affects both the equilibrium and the dynami¬ 
cal properties of the system. The spin-bath interaction 
is fully characterized by the spectral function J{uj) = 
TT^k ^k^i^ — cj/c), where we assume uj^ = Here, 

Vs represents the velocity of the sound modes of a one¬ 
dimensional Bose-Einstein condensate or a long trans¬ 
mission line. Hereafter, we shall focus on the case of 
ohmic dissipation at zero temperature, where the spec¬ 
tral function reads J{uj) = 27raujex.p Here, ujc is 

a high energy cutoff and the dimensionless parameter a 
quantifies the strength of the interaction with the bath. 
These parameters can be derived microscopically for an 
ultra-cold atom settin^^^t^. 

The bath induces both a renormalization of the tun¬ 
neling element A, and a strong Ising-type ferromagnetic 
interaction between the spins j and p, which is 

mediated by an exchange of bosonic excitations at low 


wave vectorJ^. This interaction is reminiscent of the 
Ruderman-Kittel-Kasuya-Yosida interaction for Kondo 
lattice^^. The bosonic induced-coupling has been ob¬ 
served in light-matter system^^^^^, for example. This 
interaction can be exemplified by applying an exact uni¬ 
tary transformation H = V~^HV on the Hamiltonian 


(1), with V = 

Tne transformed Hamiltonian indeed reads: 




^-k 


M 


i=i 

+ ^u)kblb, 


af + (Tj e ^ 


U-P^j^r 




( 2 ) 


where Qj = iY,k Note that = 

i+E-.i explicitly denotes the renormalized Ising cou¬ 
pling between the spins j and p, with 

\j-p\ 2 1 I ' ' 

I ^2 

The excitation of the spin j comes with a simultaneous 
polarization of the neighboring bath into a coherent state 
l^j) = |0)? resulting in a renormalization of the tun¬ 

neling element. This argument can be made rigorous 
by an adiabatic renormalization procedure, developed in 
Refs. [2] and [3l In the regime A/ujc 1, one can in¬ 
deed assume that the high frequency modes of the bath 
(above a given frequency uJi{A) corresponding to several 
units of A) adjust instantaneously to the value of the 
spin. The tunneling element is then dressed by the bath, 
and is renormalized to A < A. This procedure can be 
iterated and converges in the ohmic case and for a < 1, 
to a renormalized value of the bare tunneling element 
A to A^ = A{A/uJc)^/^^~^\ This result sheds light on 
the mechanism at the origin of the d issipativ e quantum 
phase transition induced by the batfP^^EH^ strong 
coupling, the bath entirely polarizes the spins, by anal¬ 
ogy to a ferromagnetic phase. For one spin, the quantum 
phase transition belongs to the Kosterlitz-Thouless class, 
where the order parameter at equilibrium (cr|) exhibits 
a jumpP. In this case, the critical value ac of the cou¬ 
pling is Q^c = 1. The universality class of the transition 
is unchanged when the number M of spins is increased 
(and remains finite]!^, and the a ssociated critical value 
ac decreases with (finite) due to the strong fer¬ 

romagnetic interaction between the spins induced by the 
bath. The case M = 2 was systematically studied in Ref. 

m 

The rest of the paper is organized as follows. In Sec. H, 
we summarize the general methodology used to compute 
the spin dynamics, in the case of one single spin coupled 
to bosonic degrees of freedom. These groundings will al¬ 
low us to expose the extension of the methology to the 
case of two spins (M = 2) in Sec. HI. We will then in¬ 
vestigate the quantum phase transition displayed by the 
two-spin system and present several results concerning 
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the spin dynamics both in the unpolarized and in the po¬ 
larized phase. We find that the spin dynamics in the po¬ 
larized phase exhibits an universal behaviour, in the sense 
that it becomes independent of the coupling strength a. 
Then we study the effect of the bath on the synchroniza¬ 
tion properties of the two spins in relation with spin-spin 
correlation functions. We also present Landau-Zener- 
Stueckelberg-Majorana interferometr}^^®^ protocols us¬ 
ing the interaction mediated by the bath. In Sec. IV, we 
extend the methodology to the case of an infinite array 
(M ^ oo) at a mean-field level. We present results con¬ 
cerning the dynamics as well as Landau-Zener sweeps. In 
this case, we apply a Kibble-Zurek type argument to ac¬ 
count for the mean-field dynamics. Finally, Appendices 
will be devoted to some mathematical derivations. 


II. METHODOLOGY FOR SPIN DYNAMICS 

In this Section, we re-derive the real-time spin dynam¬ 
ics in the case of M = 1 spin and introduce the nota¬ 
tions that will be used in the next sections. All the de¬ 
velopments are based on different steps related to Refs. 

andl69l which will be exposed in detail below. 

A. Feynman-Vernon influence functional 

The original reference for this technique introduced by 
Feynman and Vernon is Ref. |69l 

To compute the dynamics of the spin in contact with 
the bosonic environment, we focus on the different ele¬ 
ments of the spin reduced density matrix. Let {|cr)} = 
{\^z)A~z)} be a basis of the Hilbert spin state and 
be a basis of the bath Hilbert space e^. The total 
density matrix of the system is denoted by p, and ps is 
the spin reduced density matrix. More precisely, ps is the 
partial trace of the total density matrix over the bosonic 


degrees of freedom. The evolution of the total density 
matrix can be expressed with the unitary time-evolution 
operator of the whole system U. At a given time t, the 
elements of the spin reduced density matrix read 

{af\psit)\a'j) = '^{Un,Crf\U{t)p{to)U\t)\Un,CT'j). (4) 

n 

We have |crj), |crj) G { 1 + 2)5 \ ~z)}- Next, we express the 
propagators thanks to a path-integral description, but we 
need another hypothesis in order to go further in the cal¬ 
culations: we assume that spin and bath are uncoupled 
at the initial time to when they are brought into con¬ 
tact, so that the total density matrix can be factorized, 
p{to) = Pb(^o) ^ Ps'(^o)- For fbe remaining of the article, 
we will assume such factorising initial conditions, but the 
Feynman-Vernon influence functionnal approach can be 
generalized for a general initial condition, as shown in 
Refs. [3] and uni The initial state of the bath will always 
be a thermal state at inverse temperature f3. We start 
with the spin initially in the state 1+;^) so that 

/05(to) = |+.)(+.|=(j 0). ( 5 ) 

The time-evolution of the spin reduced density matrix 
can be then re-expressed as, 

{<^f\Ps{t)\(Tf) = J DaDa'A[(T]A*[a']T[^^^,]. (6) 

The integration runs over all spin paths a and a' such 
that |cr(to)) = \cr'{to)) = |+;,), \a{t)) = \af) and \cr'{t)) = 
jcTj). The term A[cr] denotes the amplitude to follow one 
given spin path a in the sole presence of the transverse 
field term in Eq. (1). The effect of the environment is 
fully contained in the so-called Feynm an-Vernon influ¬ 
ence functional J^[a,a'] which reads^^^ 


+’[(7, a'] = exp I - [ ds [ ds' 

I ^ jto Jto 


iT („ „o<^is)-(^'{s)cr{s')+cr'{s') , a{s) - a' (s) a{s') - a' (s') 

-iLi{s - s )-----+ L2{s - s )----- 


( 7 ) 


where a spin path jumps back and forth between the two 
values a{s) = ±1. The functions Li and L 2 read 


POO 

Li{t) = / 

Jo 

POO 

L2it) = / 

Jo 


dujJ{uj) sinci;t, 
dujJ{(jj) cos ujt coth 


puj 


(8) 


For an ohmic bath in the zero-temperature limit {3 
+ 00 , the functions Li and L 2 explicitly read. 


Li{t) = ATiauj‘1 


UJct 

(1 + uoit‘^y 


L 2 (t) = 2 t: aul ^^"^212)2 ■ (9) 

A derivation of Eq. 0 is done in the Appendix A. 

Erom Eq. 0 , we see that the bosonic environment 
couples the symmetric and anti-symmetric spin paths 
p{t) = l/2[cr(t) + cr'(t)] and = l/2[cr(t) — cr'(t)] 

at different times. These spin variables take values in 
{ —1,0,+1} and are the equivalent of the classical and 
quantum variables in the Schwinger-Keldysh representa¬ 
tion. We have then integrated out the bosonic degrees 
of freedom, which no longer appear in the expression of 
the spin dynamics, but the prize to pay is the introduc- 
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FIG. 1. (Color online) Spin states. 


tion of a spin-spin interaction term which is not local in 
time. This long range interaction in time is remini scent 
of the quantum Ising model with long range forceJI^^ 
in 1/r^. Dealing with such terms is difficult at a general 
level. The spin dynamics at a given time t depends on 
its state at previous times s <t\ the dynamics is said to 
be non-Markovian. 



FIG. 2. (Color online) Spin path- ri(t) — ~ ^j) 

red; ^{t) — ^j^{t — tj) in dashed blue. 


The prime in in Eq. (11) indicates that the initial 

and final sojourn states are nxed according to the initial 
and final conditions. More precisely we have Tq = T 2 n = 
1. The influence functional reads: 


B. “Blips” and “Sojourns” 

The next step is the rewriting of the spin path in the 
language of “Blips” and “Sojourns”, following the work 
of Ref. 121 

The double path integral in Eq. can be viewed as 
one single path that visits the four states A (for which 77 = 
1 and ^ = 0 ), B (for which 77 = 0 and ^ = 1), C (for which 
77 = 0 and ^ = —1) and D (for which 77 = —1 and ^ = 0). 
States A and D correspond to the diagonal elements of 
the density matrix (also named ‘sojourn’ states) whereas 
B and C correspond to the off-diagonal ones (also called 
‘blip’ states)!^. The four states are depicted in Eig. 

As stated previously, the spin is initially in the state 
|-h;s), so that the double spin path is initially constrained 
in the diagonal state A, which can be seen as the ele¬ 
ment top left element of the spin density matrix. We will 
first focus on the computation of the upper left diagonal 
element of the density matrix, describing the probability 

Po{t) = {+z\ps{t)\-\-z) = (1 + (^^(0))/2, (10) 

to find back the system in the state \-\-z) at time t. We 
consider then spin paths that end in the sojourn state A. 
Such a path makes 2 n transitions along the way at times 
ti, i G {l,2,..,2n} with to < ti < ^2 < ... < t 2 n- We 
can write this spin path as ^(t) = ^^0 

= 5 ^ 7=0 ~ ^ 3 ) where the variables Ei and 

take values in { — 1 , 1 }. Such a path is visualised in Eig. ^ 
The variables E (in blue) describe the blip parts, and the 
variables T (in red) on the other hand characterize the 
sojourn parts. 

After the introduction of these variables, po can be 
expressed as a series in A^, as shown in Refs. |2]and[3]: 


^n = QlQ 2 , ( 12 ) 


with 


Qi = exp 


. 2n —1 2n 

— ^j^kQl{tj ~ tk) 

k=0 j=k-\-l 


Q 2 = exp 


^ 2n—1 2n 

— ^j^kQ2{tj ~ tk) 

k=l j=k-\-l 


(13) 

(14) 


The functions Qi and ( 52 , which describe the feedbacks 
of the dissipative environment, are directly obtained from 
the spectral function J{uj) (they are second integrals of 
the Li and L 2 functions). At zero temperature, we have: 


Qi{t) = 

/ auj —^sincjt. 

(15) 


Jo ^ 


Q2{t) = 

[ (1 — cosujt ). 

(16) 


Jo 



Eor a ohmic spectral density at zero temperature, we have 

Qi{t) = 27rat8iii~^{ujct)^ (17) 

Q2{t) = 7r(alog(l + (18) 

Erom Eq. ( p^ and Eq. ( [l^ , we see that the term Qi 
couples the blips to all the previous sojourns, while Q 2 
couples the blips to all the previous blips (including self¬ 
interaction) . A derivation of these expressions is provided 
in the Appendix B. 


C. Stochastic decoupling 


Po{t) = E ( V ) / / dti y] J^n- 

' n V / Jto Jto {Ej},{rj}' 

( 11 ) 


rt2 


n=0 


At this point, the main difficulty is to treat the long 
range correlation in time induced by the bath in the 
quantum limit. The Non Interacting Blip Approximation 
(NIBA) greatly simplifies the problem and permits to 
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compute some dynamical quantities, but does not allow 
to investigate the strong coupling regime or the treatment 
of driving term^^. To decouple the spin-spin interac- 
tion(s) in time, we use Hubbard-Stratonovitch var iables 
following previous works from collaborators and 
Some efforts in this direction were also done in Refs. M- 
[34l This stochastic unravelling of the influence functional 
will allow us to write the dynamics of the spin-reduced 
density matrix as a solution of a stochastic differential 
equation. Let h and k be two complex gaussian random 
fields which verif^EB 


h{t)h{s) =-Q 2 {t - s) + /i, 

TT 


k{t)k{s) = I 2 , 


h{t)k{s) =-Qi{t - s)0{t - s) -1- Is. 

TT 


(19) 

( 20 ) 
( 21 ) 


The over line denotes statistical average, 0{.) is the Heav¬ 
iside step function and /i, I 2 and Is are arbi trary co m- 
plex constants. Making use of the identity exp(X) = 
exp(X2/2) , Eqs. (12), (13) and (14) can then be reex¬ 
pressed as: 




2n 

1=1 


+ j_i]. 


( 22 ) 


The complex constants Ip do not contribute to the aver¬ 
age because = 0- This step was 

done in Refs. [29] and IMI with the introduction of one 
stochastic field (which is valid in a certain limit, as we 
will see later), and with two fields in Ref. (ST] The sum¬ 
mation over blips and sojourn variables {Ej} and {Tj} 
can be incorporated by considering a product of matrices 
of the form 


(23) 


{A, B,C,D}. This rewriting was originally introduced 
in Ref. Then, we get 



/ 

0 

_^h+k 

0 


e' 

h-k Q 

0 

_ k 

—e 

— h—k Q 

0 

^—h+k 


V 

1 

1 

1 

0 

^h—k 

0 

the four 

dimensional 

vector 

space 


po{t) = E (V) / n 

„=0 V / Jto Jto 


Voitj). (24) 


We remark that Eq. (24) has the form of a time- 


ordered exponential, averaged over stochastic variables, 
so that we finally have: 


poit) = 


(25) 


where (4>j| = (e ^h 2 n) q ig solution of 

the Stochastic Schrodinger Equation (SSE), 


idt\^) = Vo{m) (26) 

with initial condition |4>^) = 0, 0, 0)^. 

The vector |4>(t)) represents the double spin state 
which characterizes the spin density matrix. The vectors 
|4>^) and |4>j) are related to the initial and final condi¬ 
tions of the paths. As spin paths start and end in the 
sojourn state A, only the first component of these vec¬ 
tors is non-zero. The choice of the phases is iinked to 
the asymmetry between bhps and sojourns (see Eq. (13) 
and Eq. The contribution from the first sojourn 

is encoded in |4>i), and we artificiaiiy suppress the con¬ 
tribution of the iast sojourn via |4>/). This finai vector 
depends on an intermediate time, but we can notice that 
repiacing 0, 0) by 0, 0, 0) does not add 

any contribution on average. The numerical procedure 
requires a large number of realizations of the fields h and 
k. Eor each realization, we solve the stochastic equation 
and (cr^(t)) is obtained by averaging over the results of all 
the realizations. In general we use Eourier series decom¬ 
position for the sampling of the fields h and k. Details 
about the sampling can be found in the Appendix C. 

This framework refers to as the Stochastic Schrodinger 
Equation (SSE) method. 

It is also possible to incorporate driving effects in the 
framework of the SSE method in an exact manner. In 
the present article, we will consider specifically a driv¬ 
ing term acting on the spin, which reads e{t)a^. Erom 
the path integral approach and the blip-sojourn decom¬ 
position, we see that this effect can be incorporated by 
adding a deterministic part to the stoch astic field h of 
the SSE. The resulting field readJ^^ESI 

h^{t) = h{t) + f dse{s). (27) 

Jto 


D. Previous results and discussion 


As can be seen in Eq. (23), the effective Hamiltonian 


for the spin density matrix is not Hermitian (in gen¬ 
eral h and k have both a real and an imaginary part). 
The complexness of h and k may give rise to numerical 
convergence problems at a general level, and accessing 
the regime of strong coupling between spin and bath re¬ 
quires a special attention on this issue. Eor the ohmic 
spin-boson model, simplifications occur in the regime 
A/cJc ^ 1, as shown in Refs. (29] and l3Ql In this regime 


the function Qi in Eq. (17) can be considered as a con¬ 
stant (tan“^(co’ct) 7r/2y7 allowing us to use only one 
stochastic field h which is purely imaginary. As presented 
in the references mentioned above, the SSE method then 
leads to a correct prediction of the dynamical behavior 
for 0 < Of < 1/2. 

The method with two stochastic fields was then used 
to compute the dynamics of the driven dissipative Rabi 
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modePIl, for which the use of complex fields were not 
problematic. We could in particular access large values 
of the light-matter coupling. In the previous subsections, 
we focused on the computation of the top-left diagonal 
element of the density matrix {+z\ps{'f^)\^z)^ with the 
initial condition ps{to) = \-\-z){+z\ given by Eq. 0. It 
is possible to either compute off-diagonal elements of the 
density matrix or consider another initial state for the 
spin in the framework of this method, by considering 
other initial and final vectors and |^/). Such 

developments are presented in the subsections B and C 
of the Sec. II of Ref. [3T| It is also possible to incorporate 
driving effects on the bosonic degrees of freedom, as 
shown in the subsection D of the Sec. II of Ref. [3ll 


A below. In this case, the spin-reduced density matrix 
has a dimension 16 and it is possible to develop the same 
formalism as in the one spin case, in an exact manner. 
The case of two spins is particularly interesting as the 
quantum phase transition from the unpolarized phase to 
the polarized phase occurs for a smaller value of oP^. 
While the quantum phase transition was not accessible 
with the SSE method in the case of one spin {ac = 1), it 
will be possible to investigate this regime for two spins 
(ac — 0.2), as shown in the subsection B. Synchroniza¬ 
tion is studied in the subsection C. We finally investi¬ 
gate Landau-Zener-Stueckelberg-Majorana protocols in 
the subsection D. 

A. Exact method for two spins 


Some authors did not express the spin paths in the 
language of blips and sojourns, but rather reached an 
effective stochastic Liouville equation for the density 
matrix, see Refs. I32ll35l This technique has notably 
been used to compute the dynamics for the M orse 
oscillatoi!^. Non-Markovian master equation d^^ * ^^ * were 
derived thanks to the same Eeynman-Vernon influence 
functional starting point. A review of the different 
path-integral methods developped to tackle the non- 
Markovian dynamics in spin-bath systems is provided in 
Ref. [71 

Next, we go further and present other applications of 
the method to the case of two spins (Sec. Ill), and the 
case of the array (Sec. IV). Several applications we will 
focus on have not been yet addressed in the literature 
using an alternative approach. 


III. TWO SPINS 

In this Section, we focus on the case of M = 2 spins. In 
this case, it is possible to reach an exact linear stochas¬ 
tic differential equation describing the dynamics of the 
spin reduced density matrix, as shown in the subsection 


Eor two spins, we will neglect the spatial separation 
between sites xi = X 2 = 0. We proceed as in the one- 
spin case and follow the steps exposed in Sec. II. As 
before, the two spins initially in the state 1+;^) so that 
Ps{to) = |+;2 5 ~\~z){~\~z +2 I- The time-evolution of a given 
element x = {aij^(j 2 j\ps{t)Wi f) of fho spin re¬ 

duced density matrix can be then re-expressed as. 


2 2 

P=1 P=1 

X exp |i J dsK [ai{s)a 2 {s) — a[{s)a 2 {s)]^ . 

(28) 

The integration runs over all spin paths ai, cr2, (j[ and 
such that \(Jp{to)) = Wpito)) = 1+;^), Wpit)) = \(Jf,p) 
and |<jp(t)) = As in the one-spin case, the terms 

of the form A[(Tp] denote the amplitude to follow one 
given spin path cFp in the sole presence of the transverse 
field term acting on the spin p. The last term of the 
right hand side of Eq. (28) comes from the Ising inter¬ 


action between the two spins. The influence functional 
^[ o ' 1 , 0 ' 2 , 0 - xi ^ 2 \ ^oads . 




[cr I ,(72 jO",( 72 } 


= e 


1 rt 7 rs 7 / s -^2 f r / f \ — a'-(s) cr j (s') + cr'-(s') . i\<^i 

~ ^ = l ) 2 2 \-L2{s — S ) — 


(s)-crps) CTj(s')- 




(29) 


r 


The additional term Q in Eq. (29) reads : 


Q[(7l,a2,(T[,<T2] = e 


if /; d* V ' 


of the right hand side of Eq. (28)). In the following we 
gather these two contributions in a functionnal Q which 
reads 

g[<Ti,<T2,<T[,a!2]= (31) 

with Kr = K aujc the renormalized Ising interaction. 

The paths introduced in Eq. (28) can be viewed as one 
single path that visits the sixteen states corresponding to 
tween the spins. The term above is indeed similar to the the matrix elements of the spin-reduced density-matrix, 

one coming from the direct Ising interaction K (last term We will note £ = {AA, AB, AC, AD, BA, BB, BC, BD, 


(30) 


with = 2 jir J{uj)/uj = AauJc- We recover in Eq. (30) 
that the bath renormalizes the direct Ising interaction oe- 
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FIG. 3. (Color online) Spin path for the dimer problem- The 
upper part shows the spin path in terms of blips and sojourns 
for the first spin, while the lower part shows the spin path 
of the second spin. — tj) in red; = 

~ tj) in dashed blue. The system starts in the 
state AA, jumps to the state AB at si = , then to the state 

CB at S 2 = t\. It finally ends in the state AA at t. 


CA, CB, CC, CD, DA, DB, DC, DD} the set of these 
states - the states A, B, C and D have been defined for 
one spin in Sec. II. The four states AA, AD, DA and DD 
correspond to the diagonal elements of the canonical den¬ 
sity matrix, while other states correspond to off-diagonal 
elements. As before, we consider the case where the spin 
subsystem starts in the state \-\-z^+z)^ and we intend to 
compute the probability 

Pl{t) = {-\-z^-^z\ps{t)\-\-z,-\-z)^ (32) 

to come back in the same state \-\-z^-\-z) at time 
t. Then, both the first and the second spin path 
make an even number of transitions along the way 
at times j G {1,2, ..,2np} for p G {1,2} such 
that to < 1 1 < 1 2 < ... < t 2 np < t. We can write 
these spin paths as ^^(t) = ~ ^^) and 

V^i'^) — ^^^{t — t^) where the variables and 

take values in { — 1,1}. Such a path can be visualized in 
Fig-i as a couple of one-spin paths. 

The probability pi (t) is given by a series in : 


pi{t) = 


.no V / ^ So J So 


ni,n2 


(33) 

where N = nin 2 and {^o, si,..., 52(ni+n2)} 

dered reunion of the two sequences {tj} and {tj}. The 


summation over rii and n 2 goes from 0 to infinity. The 
prime in {Tj}' in Eq. (33) indicates that the initial and 


-*-0 


= rl = 


Y'l 

2ni 


final states are fixed according to 
T2^2 — 1- The influence functional can be written ex¬ 
plicit ely in terms of Ej and Tj variables: 


•Fni,„2 = (II 1-^2 j ^[^1’^2,0-;, 0 - 2 ], (34) 


with 


Q j = exp 


Q 2 = exp 


Ail = exp 


Alj = exp 


2np —1 2np 

I Y. T. - <E) 

k=0 j=k-\-l 
2np — l 2np 

- E E 


7T 


k=l j=k-\-l 


2np — l 




j ^ k 


2nw—1 


^E E 




(35) 

(36) 

, (37) 

• (38) 


In Eqs. (37) and (38), p = 2 if p = 1 and p = 1 
if p = 2. The terms Aif and Ai^ account for retarded 
interactions between the two spins, mediated by the bath. 
Their expression in terms of blip and sojourn variables is 
very similar to the ones of Qj and Qj principle 

of their derivation is the same as in the case of one spin 
(see Appendix B). The situation differs however slightly 
since the blip variables corresponding to one spin and 
the sojourn variable corresponding to the other one can 
be simultaneously both non-zero. A detailled derivation 
in this particular case is provided in Appendix D. The 
(renormalized) Ising interaction (in ^[cri, cr2, cr^, cr^]) can 
be expressed in a convenient way in this description, as 
we have 


cri(s)cr2(s) - a[(s)a2(s) = 2 (s)C^(s) + r?^(s)C^(s)] . 

(39) 

As for the one-spin case, we can proceed to a stochastic 
unravelling of the influence functional, and we have 


2ni 


•Fni,n2 =n®^P + fc(Tl)T-l] 


i=l 


2n2 


X JJ exp [h(t])E] + 
i=i 

X 0[(Ti,C72,(T;,(T2]- (40) 


The fields h and k verify the correlations of Eqs. ( |19[ ), 
(20), and (21). Eq. (33) together with Eq. (40) has now 
























































the form of a time ordered product, averaged over the 
noise variables. 

The summation over the variables {Sj} and for 

p G {1,2} can be incorporated by considering an effective 
Hamiltonian Hi{t) for the spin density matrix, acting on 
the space £. It can be written as a sum of two terms 
Hi{t) = The (renormalized) Ising interaction 

is contained in the first term l/i, while the second term 
Vi{t) accounts for tunneling events. 

is a diagonal matrix, whose elements are {Ui)- ■ = 
2 Kr{r]l^f + T]i^i), where ryf and are the value of 77 ^ 
and for the state in the position i in the set £ = {AA, 
AB, AC, AD, BA, BB, BC, BD, CA, CB, CC, CD, DA, 
DB, DC, DD}. We sequence {Ui)- ■ gives explicitely 
( 0 ,/c, —/c, 0 ,/c, 0 , 0 , —/c, —/c, 0 , 0 ,/c, 0 , —/c,/c, 0 ) with k = 
2Kr. 

The 16 by 16 matrix Vi{t) accounts for tunneling ele¬ 
ments and has the following form. 


Finally, the dynamics of the 16 dimensional spin re¬ 
duced density matrix is governed by an effective SSE with 
Hamiltonian Hi: 

pi{t) = (43) 

where (i>/| = {e-2*(*2iv),o,0,0,0,0,0,0,0,0,0,0,0,0,0,0) 
and |i>) is the solution of the stochastic Schrddinger 
equation 

idt\<!>) = (44) 

with initial condition 

l^i) = , 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)^. (45) 

Similarly to the one-spin case, simplifications occur in 
the scaling regime, as shown in Appendix E. In this Ap¬ 
pendix, we also investigate numerical convergence issues, 
as well as other initial and final conditions, which lead to 
a different choice for the vectors and |^/). 


U(t) = f 


/ W Db^a Dq^a (0) 
Da^b W (0) Dd^b 
Da^c (0) W Db^c 

V (0) Db^d Dc^b W 


. (41) 


Each term of this matrix corresponds to a transition from 
one state in £ to another, induced by one spin-flip. It is 
written in Eq. (41) in a block structure. Each block is a 


4 by 4 matrix that can be given a physical interpretation. 
The diagonal matrices correspond to flips of the second 
spin, the first one left unchanged. As a result the matrix 
W(t) has the same structure as in the one-spin case. 


W{t) 


0 

eh-k 

_ ^—h—k 

0 


^—h+k 

0 

0 


_^ 

0 

0 

eh-k 


0 \ 

_g/l + /c 

^—h+k 

0 / 


(42) 


All the elements of the 4 by 4 matrices on the diagonal 
running from the lower left to the upper right are 
zero, because the corresponding states are not coupled 
by one single spin-flip. The eight matrices Db^a, 
Dc^a, Da^b, Db^b, Da^c, Db^c, Db^b and 
Dq^b describe spin flips of the first spin (the precise 
transition corresponds to the subscript), the second one 
left unchanged. They read respectively x /4, 

-e^+^ X /4, X /4, -e^+^ x h, x 

e~h+k X _^-h-k X and x I/i (/4 is the 

identity). Let us exemplify such transitions thanks 
to the path of Eig. The first transition at si = 
corresponds to the transition AA^AB. Its amplitude is 
given by the term of the first column and the second raw 
of the top left matrix W. The next transition at S 2 = t\ 
corresponds to the transition AB^CB. Its amplitude is 
given by the term of the second column and the second 
raw of the matrix Da^c- 


B. Nonequilibrium dynamics and quantum phase 
transition in the dimer model 

Here, we apply the SSE methodology in order to tackle 
the non-equilibrium spin dynamics in the presence of 
strong dissipative interactions in the case of two spins. 

We define the triplet subspace spanned by 
the three states {|T_) = \—z,—z)A'^o) = 

l/^/2[\^z,—z)-\-\—z,^z)]A'^-\-) ~ \-\-z,-\-z)}, while 

the singlet state is \S) = 1/\/2[\^z,-z) - \-z,^z)] 
and remains isolated in the dynamics. This problem 
is well-k nown to exhibit a dissipative quantum phase 
transitwhere the bath entirely polarizes the 
two spins either in the \T_^) or |T_) state, by analogy 
to a ferromagnetic phase. The transition line can be 
located thanks to the evolution of the entanglement 
entropy with respect to a (see Eig. 5 of Ref. EH) or 
to the evolution of the connected correlation function 
C = (erferf) “ (^i)(^ 2 ) Fig. 10 of Ref. |42j). 

With our approach, we are able to address the non¬ 
equilibrium dynamics of the system both in the unpolar¬ 
ized and in the polarized phase, and thus reproduce the 
quantum phase transition. We consider that the system 
initially starts from the state \T_^) at the time to, when 
spin and bath are brought into contact. We show in Eig.[^ 
the time evolution of P\To)^ P\t+) P\T-)^ which are the 

occupancies of the states |To) , jT^) and |T_). 

The different panels correspond to different values of 
a from a = 0.01 (top left) to a = 0.14 (bottom right). 
All these values corrrespond to the unpolarized phase in 
the range of used parameters {K = 0 and ujc = 100). 

We first note in Eig. |^a progressive suppression of the 
Rabi oscillations between the two states |Ta) and |T_) 
when increasing the parameter a. This behavior is simi¬ 
lar to the one observed in the case of the single spin-boson 
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FIG. 4. (Color online) Dynamics of the dimer model in the 
unpolarized phase: the dashed yellow line represents P|t+), 
the full red line represents P\t_) and the dotted green line 
represents P\To)- From the top left to the bottom right, we 
have a = 0.01, a = 0.02, a = 0.04, a = 0.06, a = 0.08, 
a = 0.1, a = 0.12 and , a = 0.14. The system starts in the 
state |T+) for all the plots. We have taken uod^ — 100 and 
K — ^ for all plots. 




FIG. 5. (Color online) Left panel: evolution of ln(p|T^)) at 
a — 0.14, for A/aJc = 0.01 (yellow line-top) and /S./ujc — 
0.05 (blue line-bottom), and bi-exponential fit (dashed black 
line). Right panel: Critical line with respect to /\/uc Sit K — 
0 (green dots and full green line) and comparison with the 
results obtained in Ref. [23] (TDNRG) (red triangles and red 
dashed line) and Ref. [42] (QMC) (blue squares and dotted 
blue line) 



FIG. 6. (Color online) Critical line with respect to K for 
/S./ujc — 0.01 (blue points and full blue line). Above the line, 
the system relaxes to a polarized steady-state. The red dots 
and the dotted red line show the location of the crossover 
line from coherent to incoherent behaviour for the spin oscil¬ 
lations. 


model, where the crossover from coherent oscillations to 
an incoherent dynamics occurs at ad2. At high values 
of (a, the relaxation from the initial state |T+) becomes 
slower due to the strong ferromagnetic interaction, and 
it is numerically harder to investigate the dynamics in 
the zone a > 0.1, due to the time scales involved (other 
initial states lead to an easier numerical investigation, 
allowing to determine accurately the equilibrium density 
matrix at long times). In the zone ad2 < a < ac^ we 
find a monotonic relaxation towards the equilibrium. In 
this zone, for the case of one spin, conformal field theory 
has predicted that several timescales are involved in the 
dynamics, leading to a multi-exponential deca}!^ (which 
has not been seen in NRGP^. A bi-exponential decay 
was found in this case thanks to a multilayer multicon¬ 
figuration time-dependent Hartree method^. Here, for 
two spins and at small to intermediate times, we obtain 
results which are also consistent with a bi-exponential 
relaxation, as shown on the left panel of Fig. Other 
studies have predicted more complicated forms for the 
relaxation, without any pure exponential decay (see for 
example the results of Ref. [78] obtained with renormal¬ 
ization group methods). 

We are then able to locate the phase transition 
from the divergence of the associated time scale. The 
transition line is shown on the right panel of Fig. 
together with the previous results obtained with a time 
dependent Numerical Renormalization Group (TDNRG) 
method^, or with a Quantum Monte-Carlo (QMC) 
method^. This plot corresponds to a vanishing direct 
Ising interaction K = 0, and different values of ujc- 
The phase diagram of the system with respect to the 
parameter K is shown in Fig. The full blue line 
shows the phase transition line between the polarized 
and the unpolarized phase, while the dotted red line 
shows the crossover line from coherent to incoherent 
Rabi oscillations in the dynamicJ^. 























































10 



A' (t-to) 


FIG. 7. (Color online) Universal dynamics of the dimer model 
in the polarized phase. The system starts in the no neq uilib- 
rium state described by the density matrix of Eq. (46), and 
relax towards a statistic superposition of |T+) and \T-). The 
parameters are a = 0.2, c^c/A = 100 (red points); a = 0.25, 
cJc/A = 50 (right pointing green triangles); a = 0.22, 
UcIA = 80 (left pointing blue triangles); a = 0.3, uJcj^ — 20 
(black squares). Taking K ^ 0 gives the same exponential 
relaxation. 


Next, we show results concerning the dynamics in the 
polarized phase {a > o^c), corresponding to a quantum 
quench across the critical line, from a = 0 to a > ac. 
Some theoretical studies have focused on this question 
in spin^^^M^ or bosonic system^^^^^. For example, at 
K = 0 and a = 0, the initial state of the system is given 

by iv-) = I-.) ® I-.) = 1/2(|T+) + |T_)) - l/V2\To). 

The associated spin density matrix is 


Ps{to) 


I 1 
-1 
-1 
V 1 


-1 -1 

1 1 

1 1 

-1 -1 


1 \ 
-1 
-1 
1 / 


(46) 


After a sudden change of the parameter a, the system 
is in a nonequilibrium state. We compute the spin dy¬ 
namics for different values oi a > ac and for different 
values of A/ujc- We find numerically that the system 
evolves towards the final density matrix 


lim ps{t) = I 

t—yoo 2t 


/ 1 0 0 0 \ 

I 0 0 0 0 I 
0 0 0 0 ’ 
V 0 0 0 1 / 


(47) 


corresponding to a statistical superposition of the states 
|T+) and |T_) (up to an error of around 10“^). We find 
moreover that the spin dynamics is universal in the po¬ 
larized phase, in the sense that it does not depend on a 
and K. More precisely, we find that 


Pi+-)W 


P\-+){t) =poexp 


- to) 
UJc 


(48) 



FIG. 8. (Color online) Equilibrium value of (crfcr|) as a func¬ 
tion of a for A/UJc — 0.01 (green circles), A/ujc — 0.05 (yellow 
right-pointing triangles), A/ujc — 0.1 (red left-pointing trian¬ 
gles) and A/ UJc — 0.2 (blue squares). We have A = 0. The 
lines correspond to the value predicted by a toy-model of two 
interacting spins with tunneling element A^ and Ising interac¬ 
tion Kr obtained thanks to a variational procedure. Param¬ 
eters are A/ujc = 0.01 (full green line), A/ujc — 0.05 (yellow 
dashed line), A/ujc — 0.1 (red dotted line) and A/ujc — 0.2 
(blue mixed line). The inset shows the evolution of 
with a for A/ujc — 0.2. The sign of this quantity changes 
when increasing a. 


as shown in Fig. for a quench from a = 0 to a > ac. 

_)(t) (p|_is the probability to find the system 

in the state \^z^—z) {\—z^-^z)) at time t, given by the 
diagonal term of the density matrix [ps ]22 ([pslss)- This 
simple form of the damping, and its independance with 
respect to A or a, can be accounted for by a very fast 
relaxation towards the spin ground state, without the 
emission of photons. The strong bath-induced Ising 
interaction and the orthogonality between the polarized 
state lead to a rapid evolution independent of the other 
external parameters. 

We also remark that, in the unpolarized phase, the 
value of (erferf)eg is non-zero due to the strong ferro¬ 
magnetic interaction mediated by the bath. We com¬ 
pute this quantity as the limit of tvB [ps'(^)e^ie^ 2 ] long 
times, and plot its evolution with respect to a for differ¬ 
ent values of uJc in Fig. At very small A/uje we have 
roughly {<Ji<J 2 )eq = OiUJe/ a/ {(^0 Jc)‘^ + A^, which would be 
the equilibrium value of this quantity in a two-spins Ising 
model governed by the Hamiltonian 

= ^ « + <^ 2 ) - Krcr^i^l, ( 49 ) 

where A^ = A(A/co’c)^/^^~^^ is the renormalized tun¬ 
neling element obtained by an adiabatic renormalization 
procedur^^ (see the Introduction). 
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There are notable deviations with respect to this toy- 
model, especially when A/cJc becomes larger {A/uJc > 
0.02). In this case, the adiabatic renormalization proce¬ 
dure is no longer valid, as the bath and spin degrees of 
freedom evolution time scales are not well separated. The 
assumption of fully polarized bath states associated to 
one given spin polarization no longer holds and we need 
to refine the analysis, for example by using a variational 
technique on the ground state wavefunction following the 
ideas of Refs. im and HSl We write the Hamiltonian 
of the system in a displaced oscillator basis defined by 
the four states {|-B++) (g) 1 +^, +z), \Bo) ® \+z, -z), \Bo) ® 
\-z,+z), \B —) 0 \-z,-z)}, with 


B++ = n exp 
k 

B _= 

k 


^k 


{bl - 

) 


iHn 


LUk 


\Bo). 


(50) 

(51) 


where \Bo) is the ground state of the bosonic bath taken 
in isolation at zero temperature, fk are variational pa¬ 
rameters with fk 7 ^ A/c at a general level. With this 
ansatz we do not specify the amplitude with which a 
given mode is displaced ab initio^ but these coefficients 
are found by minimizing the free energy of the total sys¬ 
tem. The displacement from the equilibrium position of 
a given oscillator may then depend on other parameters. 
Following Ref. [41] we find self-consistent equations for 
the bath-induced Ising interaction and the renormal¬ 
ized tunneling element A^, 



FIG. 9. (Color online) Synchronization phase diagram in the 
case of direct Ising interaction Kr — K. Region I (in white) 
corresponds to the unsynchronized regime : {alal) vanishes 
periodically. The yellow star shows the point for which we 
compare direct and bath-induced interaction (see text). Re¬ 
gion II (in blue) corresponds to the synchronized regime : 
(crfcr|) > 0 at all times. 


and the dissipative quantum phase transition. From the 
analytical point of view, we also note some efforts with 
multi-polaron approache^^. As seen in Fig. the main 
effect at large cJc/A is to induce a large ferromagnetic 
interaction. We will use this feature below in the syn¬ 
chronization and LZ interferometry phenomena. 


= Aexp [-a B 

I Jo w 

pOO 

Kr=a duG{u)[2 - 

Jo 


G{w) = 




(52) 

(53) 

(54) 


We plot the corresponding evolution of (crfcr 2 )eg = 

Kr!\J{KrY + A^ with respect to a for different values 
of ujc in Fig. We find a good agreement with the ex¬ 
act results given by the SSE method as long as A/cJc 
remains small {A/uJc < 1). We notably recover a change 
of the concavity of (erferf)eg with respect to er, as shown 
in the inset of Fig. where we plot the evolution of 
the second derivative of (erferf)eg for A/uc = 0.2. This 
feature cannot be recovered by the adiabatic renormal¬ 
ization procedure, but we see that this effect is far more 
pronounced in the results of the SSE than in the varia¬ 
tional treatment. The dynamical adjustment of both the 
bath and spin degrees of freedom can thus explain some 
features of the results obtained numerically, especially at 
small A/uje < 0.1 but this variational approach fail at 
quantitatively describing the regime of strong coupling 


C. Synchronization 

Synchronization phenomena occur spontaneously in a 
wide range of physical system^^. Here we quantitatively 
describe synchronization mechanisms between two spins 
1/2 starting from the polarized state without 

drive. In this two-spin problem coupled to a ohmic bath, 
some results were also obtained using the NRCf^. A 
comparison between classical and quantum regimes for 
this kind of problems without dissipation was recently 
done in Ref. |88l 

We consider the dynamics of two interacting spins with 
different bare oscillation frequencies Ai and A 2 (with 
Ai > A 2 > 0), starting from the same initial state. 
We quantify the synchronization due to the interaction, 
thanks to spin-spin correlations in time. We will com¬ 
pare the case of direct versus bath-induced interaction. 
We denote by Kr the effective strentgh of the interaction 
between the spins. In the case of a coupling through the 
bath we identify Kr = oloOc while we have Kr = K in the 
case of a direct Ising interaction. Some efforts were done 
to study this effect in Ref. |23l 

Let us first consider the case of direct Ising interaction 
K. A quantitative description of this type of synchro¬ 
nization can be done by studying the time-evolution of 
(erferf). The system starts in the state \-\-zi-\-z)i so that 
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FIG. 10. (Color online) Panels a and b: time evolution of 
(cTicrl) for a direct Ising interaction interaction Kr = K 
(panel a) and for a bath-induced interaction Kr = aujc (panel 
b). Panels c and d: time evolution of {af) and (cr|) for a direct 
Ising-like interaction Kr = K (panel c) and a bath-induced 
interaction Kr = auJc (pannel d). We have KrjAi = 0.4, 
A 2 /A 1 = 0.1 and ujc = 20Ai. 


(crfcr 2 )(to) = 1 at the initial time. We define the syn¬ 
chronized regime as the region in the parameters space 
for which (crfaf) stays positive at all times. We show in 
Fig the synchronization phase diagram with respect to 
A 2 /A 1 and iF^/Ai = KjAi. In the region I (in white), 
the two spins are not synchronized and the correlation 
function (erferf) changes sign periodically. In the other 
region (region II in blue in Fig. {<^ 1 ^ 2 ) always stays 
positive. For K^jA\ > 1/2 the Ising interaction domi¬ 
nates and the dynamics is synchronized for all values of 
A 2 . When A 2 approaches Ai, the two spins have com¬ 
parable oscillating frequencies and the synchronization is 
then easier. 

The dissipative case, for which the interaction orig¬ 
inates from the interaction with the bath, shows a 
similar phase diagram. There are however notable 
differences in the unsynchronized regime close to the 
transition line. In this region, the interaction with 
the bath leads to an effective synchronization after a 
short time unsynchronized dynamics. To exemplify this 
effect, we focus on the spin dynamics at iF^/Ai = 0.4 
and A 2 /A 1 = 0.1 in both cases. These parameters 
correspond to the yellow star in Fig. ^ The evolution of 
(cr|) and (erferf) is shown in Fig. 10 in both cases. We 
remark that in the case of direct T^ng coupling (panel 

a) , there is no synchronization transition as (erferf) 
changes sign periodically. By contrast, we remark that 
(erferf) only vanishes a finite number of times (see panel 

b) . After this short time behaviour, the system enters a 
synchronized regime for which (erferf) no longer vanishes 
and tends to a non-zero equilibrium value corresponding 
to a polarized equilibrium state. 


This synchronization effect is the sole consequence 


of the Ising-like interaction between spins. We found 
that dissipation processes enhance synchronization, as 
they favor the evolution towards more stable polarized 
states. The cases of Markovian or Non-Markovian bath 
may lead to the comparable enhancement. We note 
recent experiments in ultra-cold atoms exemplifying 
the synchronization phenomena between bosons and 
fermion^^. 


D. Landau-Zener-Stueckelberg-Majorana 
interferometry 

In this Section, we investigate the non-equilibrium 
behavior of the dimer system under an additional linear 
driving term 6(t)/2 Y^j=i 

We focus on a single linear passage, known as Landau- 
Zener problem. It corresponds to e(t) = eo ^ v{t — to)^ 
(v > 0). We choose eo < 0 with |eo|/A ^ 1 so that 
the initial state |T+) corresponds to the ground state at 
the initial time to- LandarP^, Zenei^l, StueckelberJ^ 
and Major angP^ provided an analytical description of this 
problem in the case of an isolated two-level system sub¬ 
ject to a linear sweep (K = 0 and a = 0). The survival 
probability pi^ that the spin remains in its initial state 
after the sweep, is fully determined by the velocity of the 
sweep V, and we have piz = exp[—7rA^/2'z;]. It was shown 
in Refs. M and M that the presence of a gaussian dissi¬ 
pative bath does not affect the transition probability in 
the case of the Landau-Zener sweep for one single spin, 
as long as the coupling is along the z-direction. It is no 
longer true for two spins and the presence of the bath 
affects the final state. 

For a symmetric drive only the triplet states are cou¬ 
pled to the bath, and three levels participate to the dy¬ 
namics. The system then constitutes a SU(3) Landau- 
Zener-Stueckelberg-Majorana interferometer^. 

In Fig. El we plot the different probabilities 
P\T){t 00), for |T) e {|r_), |To), |T+)} to end up in 
the state |T) at long times after a linear sweep of velocity 
V = 2A^, as a function of iF^/A. The lines correspond 
to the case where a = 0, so that the interaction between 
the two spins is only due to the direct Ising interaction 
Kr = K. We remark that the value of p\T+){t oc) 

is not affected by the Ising interaction. p\T_){t 00) 
goes to zero when the Ising interaction increases, while 
P\To){t c>o) simultaneously increases. This can be 
easily understood by the structure of the energy levels 
for the different values of K. On the upper part of 
Fig.[TTJ we draw the energy levels of the triplet states as 
a function of e, for different values of the Ising coupling 
iF, increasing from left to right. The system always 
starts at time t = to on the lower branch at negative 
bias eo- At the velocity considered here, we go from the 
regime of independent crossings (the two spins behave 
independently when K = 0, see left drawing) to the 














13 



drive, which is of great importance for quantum infor¬ 
mation purposes. 


IV. ARRAY 



0 2 4 6 8 10 12 14 

KJA 


FIG. 11. (Color online) Top: evolution of the energy lev¬ 
els with respect to the drive e, for different values of the 
direct Ising coupling at a = 0. Main figure: evolution of 
the final transition probabilities after a linear sweep of ve¬ 
locity V = 2A^ as a function of KrjA = {K + aujc)/A. 
The lines correspond to a direct Ising interaction Kr = K 
and 0 = 0, and the markers correspond to a bath-induced 
coupling Kr = auc and A = 0. Full blue line, and blue 
squares: P\T^){t oo). Dotted red line and red trian¬ 
gles: p\T_){t oo). Dashed green line and green points: 
P\To){t oo). We take Uc = 100A. 


For greater values of M, the problem becomes rapidly 
untractable numerically, as the density matrix of the 
spin system becomes too large. We will then extend the 
method at a mean field level in the case of the array 
(M ^ oc) in the subsection A. In the subsection B, 
we investigate Landau-Zener sweeps for the array and 
interpret the results with a Kibble-Zurek type argu¬ 
ment. Recent developments linked non non-equilibrium 
physics in these lattice systems involve Matrix Product 

StatesP® 

stochastic mean-field methods also allow to 
describe non-equilibrium light-matter system^^. 


A. Mean-field approximation in the limit M ^ oc 

We proceed as in the one-spin and two-spin cases and 
follow the steps exposed in Sec. II. We start with all 
the spins initially in the state 1+;^) so that ps{to) = 
n^i \^z,j){-\-z,j\- At a given time t, the elements of 
the spin reduced density matrix read 


regime of one single crossing between |T+) and |To) 
while we increase the value of K. When K/A ^ 1, the 
lowest anticrossing can be ignored and the probability 
to end up in the state |T_) then vanishes as the first gap 
closes (see the right drawing). 

The markers in Fig. correspond to the same 
protocol for A = 0 and a not zero. As can be seen, the 
dominant effect of the bath at high ujc is to induce a 
ferromagnetic Ising-like interaction. Here however, the 
probability to end up in the state |T_) does not vanish 
when increasing the value of a. This is due to transitions 
from I To) to |T_) associated to emissions of a bosonic 
excitations after the crossing of the critical point. For 
very rapid transitions, losses become negligible and the 
fidelity is higher. 

Multiple consecutive and rapid passages may result 
in constructive or destructive interferences, depending 
on the phases acquired during the adiabatic and the 
non-adiabatic evolutionJ^, allowing to propose an 
entanglement generation protocol by tuning the external 


{<^f\ps{t)\(7j) = '^{un,crf\U{t)p{to)U\t)\u„,(7j), 

n 

(55) 

where we define the M-dimensional spin vector |cr) = 

I(Ji, (72 ,.., ctm)- The time-evolution of the spin reduced 
density matrix can be then re-expressed as, 

{(^f\ps{t)Wf) = f Do-Da' exp{i[S^ - 

(56) 

The integration runs over all M-dimensional constant 
by part spin paths cr and cr' such that |cr(to)) = 

k'(^o)) = 1+) = nf=i\+z,j), kW) = k/) and 

\a'(t)) = \a'j:). Scr denotes the free action to follow one 
given M-dimensional spin path without the environment. 
This free action contains the transverse field terms, and 
the Ising interaction terms. The effect of the environ¬ 
ment is fully contained in the influence functional T[cr,<T'] ^ 
which reads in this case: 


, F ds ff ds'yz- A 


—ijCi (s — s',Xi—Xj) 


ai(s)-a'.(s) <^j{s') + a'j{s') 


-\-C2is — s',Xi—Xj) 


cT^(s) —crps) crj{s') — a‘ 


i(-'n 


X g[a,a'] 


(57) 


where we have: 


jCi{t,x) 



C- 2 {t,x) 


1 

2 


i 2 



+ L 2 



(58) 
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in the limit M ^ oo. The spins are coupled through 
three different terms: the instantaneous direct Ising in¬ 
teraction of strength K, the instantaneous interaction 
mediated by the bath in Q, and the retarded interac¬ 
tion mediated by the bath whose expression is given 
by the first term of the right hand side of Eq. (57). 
We will treat instantaneous spin-spin interactions at a 
mean field level in the thermodynamic limit M ^ oo. 
In the limit ^ 1, where a is the lattice spac¬ 

ing, we see that the retarded interactions have no ef¬ 
fect between different spins at a mean field level, since 
we have dx£i(s, x) = dx£ 2 (s, x) = 0. In the 

following, we will then neglect the retarded interaction 
between different spins, and only conserve the retarded 
self-interaction. Finally the propagation integral can be 
factorized in a product of M individual matrix elements, 
so that it is possible to write: 


FIG. 12. (Color online) Top: Evolution of the direct Ising 
interaction which is induced by the presence of the bath be¬ 
tween two spins distant of x, as a function of xujclvs. Bottom: 
Space-time dependency of the coupling functions £i (left) and 
£2 (right). The bath induces a long-range interaction between 
spins. 


The bosonic environment couples the symmetric and 
anti-symmetric spin paths r]^{t) = l/2[(jp{t) crp{t)] and 
C^{t) = l/2[cFp{t) — crp{t)] at different times and different 
lattice sites. In Fig. we plot the space and time 
coupling functions £i (bottom left) and £2 (bottom 
right). We see that the bosons induce a long-range 
interaction between spins. The maximal effect between 
two spins separated by a distance x occurs after a time 
x/vs^ due to the finite sound velocity Vs of the excitations. 


The last term of Eq. 


reads 


Qlcr,(r'] 

with IX 
the bath 
magnetie 
K' — 

expression 


of Fig. 
where x 
i and j. 
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(59) 

= 2 / 7 r Jq J{uj)/uj. We recover that 

is responsible for an indirect ferro- 
Ising-like interaction between the spins 
l/(27r) J{uj)/uj cos[{xi — Xj)/vs]^ whose 
is given in Eq. We plot on the top panel 
the value of with respect to xbod^s^ 

Xi — Xj is the distance between the two sites 


The bath is responsible for two distinct types of 
interactions. The first one is a retarded interaction 
mediated by the bosonic excitations, which travel at the 
speed Vg. The second one is an instantaneous interaction 
iF', of which we have given a physical interpretation 
thanks to the polaronic transformation in Eq. (2). 

Dealing with the spatial extent remains difficult, and 
we will treat the array problem at a mean-field level 


{^pj\ps,p{t)\o'pj) — J DapDapAp[ap]Ap[ap] J'plo'p^o'p] 

(60) 

where ps,p denotes the density matrix of spin p. Ap[ap] 
denotes the amplitude to follow a given path for the spin 
p in the sole presence of the transverse field. We have 
Kr = K Kp The remaining term J^picTp^cTp] 

encapsulates the effect of the bosonic bath on the spin p. 


Xp[<7p-,(y']=e-xpi [ ds f ds'-Li{s - s')^P{s)riP{s') 

Jto Jto ^ 

--LXs-s')e{s)e{s')}. 

IT J 

(61) 


We will drop the p index in the following, as all the sites 
are equivalent in the mean-field description. Following 
the same steps as in Sec. II, we focus on the computation 
of P 2 {t) = (+;s|ps'(t)|-h;s)and reach the same expression 
than for po{t) (see Eq. (pT])), with 

^n = SlQ 2 Q 3 . (62) 


The expressions (13) and (14) for the expressions of Qi 
and Q 2 are still valid, and we have the additionnal term 


Qs = exp 



ds{a^{s)) 


(63) 


We then reach for P 2 {t) the same expression as the 
one obtained for po{t) in Eq. (25), with the same 
final vector and |0) solution of the SSE (26), with the 
effective Hamiltonian given by (23) provided that we 
add to the stochastic field h the field hp defined by 
hi{t) = —2iKrjl^ds{(7^{s)). We have then reached 
an auto-coherent equation, as (cr^(t)) enters in the 






















15 


expression of The numerical procedure requires 

a larger number of realizations of the field h and k 
compared to the one-spin case. For each realization, we 
solve the stochastic equation and (cr^(t)) is obtained by 
averaging over the results. The effect of (cr^(t)) in hi{t) 
is dynamically updated with the number of samplings. 

We can use our method to compute the free spin dy¬ 
namics in the limit of M ^ oo. We check that the bath 
causes a decay towards one of the two equilibrium states 
in the ferromagnetic phase as well as a renormalization 
of both the tunneling element and the Ising coupling. 
However, it does not affect the university class (second 
order with mean-field exponents for the paramagnetic- 
ferromagnetic transition) of the quantum phase transi¬ 
tion as long as the direct Ising term K is not zero. This 
behavior can be understood thanks to a thermodynamic 
analysis of the action at low wave-vectors q and low fre¬ 
quency cj, which is dominated by the contribution of the 
long range Ising interaction, as shown in Appendix F. 


B. LZ transitions : Array 


We focus now on many-body Landau-Zener sweeps for 
the array, at a mean-field level. Let us underline that 
this protocol is different from the dynamical transition 
of the quantum Ising model in transverse field with near- 
est neighbours interactions studied in the litteratur^^^^^ 
(and references therein), where the driving parameter is 
the transverse field and which can be studied elegantly 
in k space. Here, we are interested in the dynamics of 
local spin variables at a mean field level. A rigorous de¬ 
scription of the dynamics should involve all the energy 
levels of the system, and their respective avoided cross¬ 
ings. Our mean-field description greatly simplifies the 
problem and the interplay of all the levels is reduced 
to a single avoided crossing governed by the local self- 
consistent Hamiltonian, 


H — — 
Hj - 2 




+ bk^ + LOkb\bk 


(64) 


The presence of the Ising interaction or the presence 
of the bath both lead to a change in the final value of 
(cr^(t ^ oo)). The origin is the same in both cases at 
weak coupling and large ujc ^ A, as the dominant ef¬ 
fect of the bath is to induce a ferromagnetic Ising-like 
interacti on K^ . In the following we use a Kibble-Zurek 
argumenfP^lHl order to describe quantitatively this ef¬ 
fect. The single site fast Landau-Zener transition can 
indeed be described thanks to the Kibble-Zurek mecha¬ 
nism, which predicts the production of topological defects 
in non-equilibrium phase transit ions^^J^. This descrip¬ 
tion splits the dynamics into three consecutive stages: it 
is supposed to be adiabatic in the first place, then evolves 




FIG. 13. (Color online) Left: schematic interpretation of the 
Landau-Zener sweep for the array in the framework of the 
Kibble-Zurek mechanism. The line [AD) shows the evolu¬ 
tion of the bare bias field with respect to time, while the 
broken line connecting points B and C represents the effec¬ 
tive bias field. The lines are full during the adiabatic stages, 
and dashed during the frozen (non-adiabatic) period. Right: 
Fast sweep (f/A^ = 8) in the array, for different values of a 
corresponding to = 0 (red curve), Kr — 2 (green curve), 
Kr — 4 (yellow curve), Kr — 6 (blue curve), Kr — 8 (ma¬ 
genta curve) and Kr — 10 (cyan curve). We have K — and 
ujc = 100A. The dashed line is the usual theoretical predic¬ 
tion for one spin, following Landau, Zener, Stueckelberg and 
Majorana; see Sec. Ill D. Inset: the blue points show the 
values of {<J^[t oo)) with respect to Krj corresponding 
to the parameters of the main plot (interaction mediated by 
the bath Kr — oluJc [K = 0)). Green squares correspond to 
a direct Ising interaction Kr — K. The full (dashed) red line 
shows the expectation value of (cr^(t ^ oo)) with respect to 
KrjlA [KrjlAr) dcduccd from the Kibble-Zurek mechanism. 


in a non-adiabatic way near the transition point, and fi¬ 
nally becomes adiabatic again. The impossibility of the 
order parameter to follow the change applied on the sys¬ 
tem provokes this non-adiabatic stage, where the dynam¬ 
ics is said to be “frozen”. It is convenient to introduce 
the characteristic energy scal^^ 

e = A/y2 { [1 + 16t>V(7r^A^)] , (65) 


which sets the limit between adiabatic and frozen stages 
(see left pannel of Fig. 13). 


We first focus on the case where a = 0 and the di¬ 
rect Ising interaction K is not zero, so that Kr = K. 
The effective field felt by one site is the sum of the bias 
field e[t) and the Ising interaction, and will be denoted 
eeff[t). The dynamics always enters in the frozen stage 
with (cr^) I, so that we have eeff{t) = e[t) — Kr during 
the first adiabatic stage. At the end of the frozen stage, 
the spin expectation value has changed, and the effective 
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field becomes eeff{t) = e(t) — Kr{(T^{t)). This leads to a 
change of the effective speed at which the frozen zone is 
crossed through, and ultimately of the transition prob¬ 
ability. This can be seen on the left pannel of Fig. 
where we show the evolution of both the bare and the ef¬ 
fective bias fields with respect to time. We can estimate 
the renormalization of the effective speed self-consistently 
thanks to basic geometrical considerations in the trape¬ 
zoid (ABCD) of the Fig. (left panel). The effective 
crossing speed is given by 


jv) + ijVeff) 

tc{Veff) - ts 


( 66 ) 


The denominator can be simplified by writing that 
— ts = — to] + {to “ ^a) “ {ts “ 

t^). We know that {to — tA) = [e{v) ^ e{veff)]/v^ 
and \tc{'^eff) — ^d] — {ts — tA) can be expressed as 
—Kr [1 — {a^{tc,Veff))]. Next we suppose that we can 
approximate {a^{tc^Veff)) by (cr^(t ^ oo)). Altogether, 
we get 


due to the presence of the bath. 

We have studied dissipative sweep protocols in the case 
of the dimer model and for an infinite array. The cou¬ 
pling to a dissipative environment is responsible for both 
an Ising-type interaction and relaxation processes. In 
the regime where A/cJc ^ 1 the predominance of the 
bath-induced Ising like interaction on relaxation mech¬ 
anisms renders possible a quantitative prediction of the 
dynamics. In the case of two spins, it was indeed possi¬ 
ble to understand the evolution of the energy levels and 
take into account the three avoided crossings. Increas¬ 
ing the number of spins would lead to a larger number of 
level crossings and a more complex energy level structure, 
as one should take into account the side-by-side avoided 
crossings of all the energy levels (except the eventual sin¬ 
glet which remains isolated). In the case of a large num¬ 
ber of spins with long-range interactions, we recover a 
local spin-1/2 view on the dynamics which can be under¬ 
stood with a single-crossing view, by analogy with mean- 
field methods. In this case, the self-consistency equation 
comes from a Kibble-Zurek type criterion with adiabatic- 
ity considerations. 


Ml = _ e(^) + ejveff) _ 

V e{v) + e(t>e//) - 2Kr [1 - PlziVeff)] ■ 


V. CONCLUSION 


It allows US to know the variation of the effective 
speed Veff at which the transition is crossed with re¬ 
spect to the Ising interaction Kr. The spin expectation 
value (cr^(t ^ oo)) is then estimated thanks to the 
Landau-Zener formula, and its evolution with respect to 
Kr = K is shown by the red curve in the inset of the 
right part of Fig. The estimation matches well the 
results obtained numerically (green squares). 


Now we take K = 0 and a not zero so that Kr = aujc. 
We plot on the right panel of Fig. the dynamics 
obtained with the SSE. We see in the inset that, at 
small Of, the estimation of the final value of the spin 
variable thanks to Eq. (67) is correct. However, it 
breaks down when the dissipation strength is increased 
because the assumption {a^{tc^Veff)) — {cF^{t oo)) 
used to derive Veff is no longer correct. Relaxation 
processes occur after the crossing of the frozen zone 
which lower (cr^(t ^ c>o)). This can be seen on the 
behaviour of the curves obtained at large values of a 
(the cyan curve for example), where the spin expectation 
value continues to go down during a rather long time 
after the crossing. The dotted red curve takes into 
account the renormalization of the tunneling frequency 


To summarize, we have developed a stochastic ap¬ 
proach to address the spin dynamics in dissipative quan¬ 
tum spin arrays. Using complex gaussian random fields, 
we have carefully studied the applicability of the method 
in all the applications. We focused first on the quan¬ 
tum phase transition displayed by two spins in contact 
with the same ohmic bath, and studied quenched dynam¬ 
ics both in the unpolarized and in the polarized phase. 
We also investigated quantitatively bath-induced syn¬ 
chronization phenomena occurring in this system. Then 
we considered non-equilibrium sweep protocols, both in 
the case of two spins and for the quantum Ising chain 
with long-range forces. In this latter case, the dynam¬ 
ics can be understood thanks to a simple Kibble-Zurek 
type argume nt. O ur results can be tested in ultra-cold 
atom systemJ^^I^. The metho d could also be applied to 
the sub-ohmic spin-b oson model^^^^^, Jaynes-Cummings 
or Ra bi a rray j^^ * ^^^ *, for topological probl ems with Dirac 
point 4^^, and for fermionic environmentas in 
Kondo lattice^^. 
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Appendix A: Feynman-Vernon influence functional 

Here we derive the expression (7) given in the main text, using the method of Ref. 11041 In order to simplify the 
derivation, we first consider that one single bosonic mode is coupled to the spin, and we have the Hamiltonian 

H + ^{b + b^)a^ +ujb^b. (Al) 

The general case of several modes will be deduced from this simpler case at the end of this appendix. Let us call 
Hs = A12(7^ the spin-part and Hb = hujb^b -j- + 6 ^) the interaction part. From Eq. (4) of the main text and 

after the introduction of the identity both on the left and on the right of the term p(to), we get 

^ \^{Un,Crf\U{t)\Um,<Jk){Um,<yk\p{to)\Up,ak'){Up,ak'\U{t)\Um,Cr'f)y (A2) 

n,m,p,k,k' 

Next, we use the factorising initial condition p{to) = psito) ^ \-\-z){+z\ and reach 


{af\psit)\a'j) = {Um\pB (^o) \ U (t) \ Ujyi^ z ) z\U (t) \ Ujyi^ (Jj) . 

n,m,p 


(A3) 


The last two terms can be expressed thanks to a path integral. The resulting action can be divided into two parts 
Ss and Sb^ the first one resulting from the spin Hamiltonian alone and the second one resulting from the remaining 
part. Factorising the spin part, we reach the equation ( 6 ) of the main text, where we get 

F[a,a'] = tVB \^pB{to)UBW]{t)Ul[a']{t)^ , (A4) 

where UBlcf] being the time evolution operator related to where a is a classical time-dependent spin-path. In 
order to evaluate this functional we need to derive the expression of the bath evolution operator. To do so, we switch 
to the interaction picture (where V = A/ 2 (a + a^)a is the interaction term) and define 1 /bM(^) the corresponding 
time evolution operator. We have : 

ihdtUB[(T]{t) = V{t)UB[(T]{t) (A5) 


Defining X = and P = the commutation relations gives: 

^ i^-f^ ^—iujb^btp^iujb^bt 

^—iujb^btppujb^bt — p _ ^—ioob^ bt pujb^ bt 

which results in: 


V(t) = —a^(t)[(b + P) cosuut + -— 
2 i 




jt]. 


(A 6 ) 


As the evolution operator pBlPiP is unitary, we suppose that we can write it as e e ^^(t)(b-\-b'^) ^ ^ 

The Schrodinger equation gives us the expression of a, /3^ 7 : 


' P{t) = ds^a{s) cos Co’S 
< 7 (t) = ds|cr(s) sincjs 

a{t) = — ds /q ds' (1)^ a{s)a{s') cos ujs' sincjs 
Then, we have : 


F[a,a'] = e*[“'(*)-“(*)] J dX{X\pB{0) ^-m)x (A7) 

where the states \X) represent a complete set of position eigenstates. It simplifies into 

F[(t,c 7'] = e*[“' J dA(A|pB( 0 )|A + 7 (i)(A 8 ) 
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In order to evaluate the element {X\pB{to)\X + — 7 '(t)), we assume a thermal equilibrium at inverse temper¬ 

ature P for the operator Pb(^o)- 


(Xi|pb( 0 )|X 2 ) = - 


g 2 sinh /3 u;q 


[[xl+xl) cosh/3a;o- 2 X 1 X 2 ] 


Z \ 27r sinh ydcjo ^ 

Using the properties of Gaussian integrals, as well as the identity ~ tanh^dcjo/^, we get: 

F[cr, a'] =e* [^'d)-«d)] +* [/3'd)-/3d)] [7d)+7'd)] coth/3a;o/2 [(/3'(t)-/3(t))^ + ( 7 (t)- 7 '(t))^] ^ 


(A9) 


(AlO) 


Hence re-inserting the expressions of a, /3 and 7 and after trigonometric calculations and using the symmetry of the 
integrand we finally recover Eq. (7) of the main text, with Li{t) = TrA^sincjot and I/ 2 (t) = ttA^ coscjot coth;dci;o/2. 
The generalization to an infinite number of modes is straightforward. 

Appendix B: Blip-sojourn development and derivation of Eq. ( |13| ) and Eq. 

Given a path (see Fig. | 2 ]of the main text for example), we can evaluate Eq. 0 of the main text. First we evaluate 
the contribution given by Li. 


rt ns ■ nt 2 j rt 2 k + l 

ds ds'Li{s — s')^{s)r]{s') =— ^ ^jrjk / ds ds'Li{s — s') 

dto Jto ^ — Q dt2j-l dt2k 


— — ^ ijdk[Ql{t2j-l — hk) + Qlihj — ^2/c+l) — Qlihj — hk) — Qlihj-l — ^2/c+l)] 

j>k=0 

2 n 

= “ E (Bi) 

j>k=0 


with Qi the opposite of the second integral of Li with Qi(0) = 0. Then we evaluate the contribution given by L 2 . 

-I -j ^ nt2j rt 2 k 1 rt2j rs 

-/ ds ds'L 2 {s - s')C{s)C{s') = -y] / ds ds'L 2 {s - s') - 77 / ds ds'L 2 {s-s') 

Jto Jto Jt2j-1 Jt2k-1 ^ Jt2j-1 Jt2j-1 

1 "" 

= — / , ^j^k[Q2{t2j-l — t2k-l) + Q2{t2j — t2k) “ Q2{t2j “ ^2/e-l) “ Q2{t2j-1 “ ^2/c)] 

TT ' 
j>k=0 

^j^jQ2{t2j ~ ^2/c-l) 

TT ^' 

J 

^ 2n 

= / ^j^kQ2{tj ^/e): 


(B2) 


j>k=0 


with Q 2 the second integral of L 2 with Q 2 ( 0 ) = 0. We recover Eq. ( 12 ), (13) and (14) of the main text, where Qi 
contains the coupling of the blips to all the previous sojourns, and Q 2 contains the coupling of the blips to all the 
previous blips (including self-interaction). 


Appendix C: Sampling of the stochastic variables and numerical convergence 


In order to sample the variables h and k which verify the correlations of Eq. (19), ( [ 2 Q| ) and ( [ 2 T] ), we use a Fourier 
series decomposition of the functions Qi and (52- To do so, we introduce the variable r = t/tf where t/ is the final 
time of the experiment/simulation. Hence r Q 2 {rtf) and r Qi{Ttf)0{T) are defined on [—1,1]. We extend their 
definitions by making them 2-periodic functions and it is then possible to expand them in Fourier series. In particular, 
we have: 
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^2 [{Tj 90 I 9 m r I / \ j,* / \ i 1 

- — - - = y + ^ [(Pm{Tj)(f>m{n) + h.C.] 


Ql [{Tj - Tk)tf] 


2^2 

m=l 

^{Tj - '^fe) = y + £ "Y {(t>m{Tj)4>*m{Tk) + h.C.] + ^ ^ [(t>m{Tj)4>*m{Tk) - h.C.] , (Cl) 
m=l 


m=l 


where (p^n • ^ exp(im7rr), and we have for m > 1, |^m = cosm7rr|, 

|/m ~ cosm7rr|, and |/^ = 0{r) sinm7rr|. go and /o are the constant Fourier 

coefficients. Then we define h and k as 


= ^^rn{rtf) 
m=l 


(si.TO + (Wl,m + ^W2,rr») + ) (I'l.m + W2,m) 

L (^) {U3,m + iU4,m) + (^ ] {V3,m + iV4,m) 


(C2) 


~ ^ ^ 4^m I “1“ T ( I “1“ 


m=l 


^ y ('^3,m “1“ '^'^4,m) y ^ ) ("^3,771 H“ '^'^4,m) 


(C3) 


where {gi,m|? {uj^rn} and {vi m} are standard normal variables. One can check that h and k verify the correlations 
given by Eqs. (p^, (20) and of the main text. In general these fields are complex, and the presence of a non-zero 


real part may lead to an exponential slowing down of the convergence. We can check that gm < ^ foi* all m, such 
that the part of the field h coming from Q 2 is purely imaginary. On the other hand, we always have terms of the 
form ui + iu 2 when it comes to the decoupling of Qi. It is then impossible to constrain the real part of the fields 
coming from the Qi decomposition. In the numerics, we use Fast Fourier Transform in order to increase the speed 
of the numerical procedure. We could have decomposed the fields on another basis of functions, but the choice of 
Fourier series decomposition seems natural, given the form of the functions Qi and Q 2 in Eq. (15) and Eq. (16). 


Appendix D: Expressions of Alf and AI 2 

The derivation of Qf and Q 2 can be found in the Appendix B. In this case, the blip and sojourn variable cannot 
be simultaneously both non-zero. Eor and fhe situation is different as the state of the first spin does not 
constrain the state of the second one. More explicitly, for Alf for example, one of the spins may be in a blip state 
while the second one is in a sojourn state, as illustrated in Eig. In the following, we will compute the contribution 
of these particuliar blip-sojourn configurations. 

The first case (left panel) yields. 


7 — j — r^2k-\-i r^2j r^2k+i 

-/ ds ds'^P{s)gP{s')Li{s-s') = -/ ds ds'Li{s-s')^ ds ds'Li{s-s') 

^ Jt^. , JtS, ^ JtP. . Jtl, As,,. Jtl. 


TT L 


^2i-1^2/cQl(^2j-l ^2/e) +'^2i^2/eQl(^2i ^2/c) + ^2/c+lQl (^2^ ^2/c+l) * 


The second configuration gives. 


-- [ ds [ ds'^^{s)gP{s')Li{s - s') = [ ds [ ds'Li{s-s') 
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The third configuration gives, 
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2fc + l *^2^ 
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The fourth configuration gives, 
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FIG. 14. (Color online) Coupling of a blip of the spin p with a simultaneous sojourn of the spin p. There are four distinct 
configurations. 


We finally recover the expression (37) of the main text. An analog computation permits to find back Eq. (37). 


Appendix E: Scaling regime 


In the scaling regime A/cJc 1, it is possible to overcome the sign problem naturally arising in our method as shown 
in Ref. [30| Simplifications occur in Eqs. (35) and (37) as we can consider that Qi{tj — tk) = 27r(atan“^ [^c{tj ~ ^/c)] — 
TT^Qf. Then we have 


2n„ —1 




lira 


^=0 


2na 




i=i 


(El) 


for g = p or p = p. is the value of ^^{t) in the interval [tj, tj+i] and pj is the value of p^{t) in the interval [tf, 

The integer I is defined by In the case of p = g, we just have I = j — 1. 

This expression does not depend on intermediate times, but only on the path taken. As a result, there is no need 
to introduce the time-dependent field k. After having introduced the field h as in the main text, we finally recover 
Eqs. (43) and (44) of the main text, with 
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where 

a = 

exp(i7rQf), 

10.)^ = 


= (1,0,0,0,0,0,0,0,0,0,0, 

JZ) 

0 

0 

' 0 

It is 

also 

possible 

to 

com- 


pute for example the probability to arrive finally in the state 1+;^,— 2 )- This can be done by taking 
10/)^ = (0, 0, 0,1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0). Similarly, one can compute the dynamics for another initial 
state. One can consider for example an initial density matrix given by Eq. ( [4^ of the main text. This corresponds 

to 10,r = 1/4(1, -1, -1,1, -1,1,1, -1, -1,1,1, -1,1, -1, -1,1). 

We can use the scaling regime simplification exposed above, even when we do not have A/cJc ^ 1- We write 
Qi{t) = 7r‘^a + [Qi{t) — 7r‘^a] and we take into account the constant part as exposed above. The remaining part 
[Qi{t) — is then decomposed into Fourier series. 

As we use a Fourier decomposition, we choose the same discretization step in time and in frequency, and take 2^ 
points. In Fig. we show the numerical convergence concerning the dynamics of p\T+){t) for the dimer problem 
with initial conmtion |T+) (see III B of the main text), with a = 0.02, ujc = 100, K = 0, for N from 6 to 11. For 
N > 11, all the curves give the same result (superposed to the full black curve). In the regime a > ac/2, one finds 
the existence of a “sweet spot” which links the final time of the simulation and a, for a given discretization. 



FIG. 15. (Color online) Time evolution of p\T^){t) for the dimer being initially in the state |T+), for A = 6 (dashed blue line), 
A = 7 (dotted green line), A = 8 (full yellow line), A = 9 (dotted red line), A = 10 (dashed purple line), and A = 11 (full 
black line). Parameters are a = 0.02, cUc = 100, and A = 0. 
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Appendix F: Thermodynamic analysis of the action for the dissipative Ising model in transverse field 

The mean-field dynamics is not affected by the presence of the bath. This behavior can be understood thanks 
to a thermodynamic analysis of the action at low wave-vectors q and low frequency cj, which is dominated by the 
peaked contribution at g = 0 of the long range Ising interaction. Using a mapping to a classical Ising model, it 
is possible to estimate the spin-spin coupling due to the environment, by focusing on the partition function (path 
integral approach) and tracing out the environmental modes 


J D{b,b*)e-^ = 


exp 


1 r/5 r/5 foo 


e ^ I + 2n^(co’) coshco’(T — r') cos ( ) aj(r)ar{r') 


B(T — T',Xj —Xr) 


(FI) 


where the aj are the classical spin variables corresponding to the eigenvalues of the quantum operators cr|, and r is 
the imaginary time. At zero temperature, we have 


B{r — U, Xj — Xr) =Re 


27raujl 

(l + UJc\t - t'\ + j 


(F2) 


where ^ = VsjoOc^ then modifying the coupling between the spins. On the other hand, the direct Ising coupling is 
responsible for a coupling term of the form 


C{t - t', Xj - Xr) = ^S{t - t'), (F3) 

and the constant behavior in the space domain dominates in the low g, low uj expansion of the action. 

The mean-field coupling then dominates over the dissipative effects and we find back the characteristic features of 
the mean-field transition of the quantum Ising model in transverse field. This mean field behavior is valid as long as 
the direct Ising term K is not zero. 
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